{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "# EGS_inf "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:35:09.142361Z",
     "start_time": "2023-11-04T11:35:08.958850Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt \n",
    "from scipy.optimize import fsolve\n",
    "from scipy.interpolate import PchipInterpolator\n",
    "import warnings\n",
    "warnings.simplefilter('ignore', RuntimeWarning)\n",
    "plt.rcParams['font.family'] = 'Times New Roman'\n",
    "plt.rcParams['font.size'] = '16'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:10.885898Z",
     "start_time": "2023-11-04T11:41:10.871938Z"
    }
   },
   "outputs": [],
   "source": [
    "r=0.05;\n",
    "mu0=0.01;\n",
    "mu1=0.03;\n",
    "sigma=0.25  ## 0.15\n",
    "alpha=0.5;\n",
    "taue=0.4;\n",
    "taui=0.05;\n",
    "I=20;\n",
    "p=18;\n",
    "x=1; \n",
    "pi0=0.5;\n",
    "eta=1;\n",
    "t0=0;  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:12.897726Z",
     "start_time": "2023-11-04T11:41:10.890886Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate some intermediate values\n",
    "psi = mu1 - mu0\n",
    "mu = mu0 + psi * pi0\n",
    "kappa = psi * (mu1 + mu0 - sigma**2) / (2 * sigma**2)\n",
    "d=r-mu1\n",
    "beta1 = 1/2 - mu0/sigma**2 + np.sqrt((1/2 - mu0/sigma**2)**2 + 2*r/sigma**2)\n",
    "beta2 = 1/2 - mu0/sigma**2 - np.sqrt((1/2 - mu0/sigma**2)**2 + 2*r/sigma**2)\n",
    "beta3 = 1/2 - mu1/sigma**2 + np.sqrt((1/2 - mu1/sigma**2)**2 + 2*r/sigma**2)\n",
    "beta4 = 1/2 - mu1/sigma**2 - np.sqrt((1/2 - mu1/sigma**2)**2 + 2*r/sigma**2)\n",
    "if np.isclose(kappa, 0):\n",
    "    kappa = int(kappa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def XB_Default_B(mu,c):\n",
    "    beta_neg = 1/2 - mu/sigma**2 - np.sqrt((1/2 - mu/sigma**2)**2 + 2*r/sigma**2)\n",
    "    xb=(beta_neg/(beta_neg-1))*((r-mu)/r)*c\n",
    "    return xb\n",
    "def E_Equity_Mu(xi,mu,c):\n",
    "    beta_neg = 1/2 - mu/sigma**2 - np.sqrt((1/2 - mu/sigma**2)**2 + 2*r/sigma**2)\n",
    "    xb=XB_Default_B(mu,c)\n",
    "    E = (1 - taue) * (xi / (r - mu) - c / r + (c / r - xb / (r - mu)) * (xi / xb)**beta_neg)\n",
    "    return E\n",
    "def D_Debt_Mu(xi,mu,c):\n",
    "    beta_neg = 1/2 - mu/sigma**2 - np.sqrt((1/2 - mu/sigma**2)**2 + 2*r/sigma**2)\n",
    "    xb=XB_Default_B(mu,c)\n",
    "    D = (1 - taui) * c/r + ((1 - alpha) * (1 - taue) * xb / (r - mu) - (1 - taui) * c/r) * (xi / xb)**beta_neg\n",
    "    return D\n",
    "def Phi_GCost_Mu(xi,mu,c):\n",
    "    D = D_Debt_Mu(xi,mu,c)\n",
    "    phi = 1 - r * D / ((1 - taui) * c)\n",
    "    return phi\n",
    "def F_Option_Mu(xi,mu,c,x0):\n",
    "    beta_pos = 1/2 - mu/sigma**2 + np.sqrt((1/2 - mu/sigma**2)**2 + 2*r/sigma**2)\n",
    "    D = D_Debt_Mu(xi,mu,c)\n",
    "    phi = 1 - r * D / ((1 - taui) * c)\n",
    "    E = E_Equity_Mu(xi,mu,c)\n",
    "    phi = Phi_GCost_Mu(xi,mu,c)\n",
    "    F = (E - phi * (1 - taui) * c/r - (I - (1 - taui) * c/r)) * (x0 / xi)**beta_pos\n",
    "    return F\n",
    "def P_FGap_Mu(xi,mu,c):\n",
    "    phi = Phi_GCost_Mu(xi,mu,c)\n",
    "    P=(1-phi)*(1-taui)*c/r\n",
    "    return P\n",
    "def Pit_2Val(xi,tc,x0):\n",
    "    pit = 1 / (1 + (1 - pi0) * ((x0 / xi) ** (psi / sigma ** 2)) * np.exp(kappa * tc) / pi0)\n",
    "    return pit\n",
    "def Phi_GCost_Pit(xi,mu,c):\n",
    "    D0t = D_Debt_Mu(xi,mu0,c)\n",
    "    D1t = D_Debt_Mu(xi,mu1,c)\n",
    "    phit = 1 - (r * ((1 - mu) * D0t + mu * D1t)) / ((1 - taui) * c)  \n",
    "    return phit\n",
    "def P_FGap_Pit(xi,mu,c):\n",
    "    phit = Phi_GCost_Pit(xi,mu,c)\n",
    "    Pt=(1-phit)*(1-taui)*c/r\n",
    "    return Pt\n",
    "def smxy(xin, yin, nx=500):  \n",
    "    Fin = PchipInterpolator(xin, yin)\n",
    "    x_sm = np.linspace(xin[0], 1, nx)\n",
    "    y_sm = Fin(x_sm)\n",
    "    return [x_sm, y_sm]\n",
    "def successive_approx(T, u_0, tolerance=1e-6, max_iter=10000, print_step=25):\n",
    "    # Computes an approximate fixed point of a given operator T via successive approximation.\n",
    "    u = u_0\n",
    "    error = np.inf\n",
    "    k = 1\n",
    "    while error > tolerance and k <= max_iter:\n",
    "        u_new = T(u)\n",
    "        error = np.max(np.abs(u_new - u))\n",
    "        if k % print_step == 0:\n",
    "            print(f\"Completed iteration {k} with error {error}.\")\n",
    "        u = u_new\n",
    "        k += 1\n",
    "    if error <= tolerance:\n",
    "        print(f\"Terminated successfully in {k} iterations.\")\n",
    "    else:\n",
    "        print(\"Warning: hit iteration bound.\")\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.array([0.8, 0.8])\n",
    "def FB1(P0):\n",
    "    xi1, c1 = P0\n",
    "    xb1 = XB_Default_B(mu1,c1)\n",
    "    F1 = ((taue - taui) * c1/r - I) * beta3 + (beta3 - 1) * (1 - taue) * xi1/(r - mu1) + (beta3 - beta4) * ((taui - taue) * c1/r - alpha * (1 - taue) * xb1/(r - mu1)) * (xi1/xb1)**beta4\n",
    "    F2 = (1 - taui) * c1/r + ((1 - alpha) * (1 - taue) * xb1/(r - mu1) - (1 - taui) * c1/r) * (xi1/xb1)**beta4 - p\n",
    "    return [F1, F2]\n",
    "xi1, c1 = fsolve(FB1, y0)\n",
    "xb1 = XB_Default_B(mu1,c1)\n",
    "E1 = E_Equity_Mu(xi1,mu1,c1)\n",
    "D1 = D_Debt_Mu(xi1,mu1,c1)\n",
    "phi1 = Phi_GCost_Mu(xi1,mu1,c1)\n",
    "Fg = F_Option_Mu(xi1,mu1,c1,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.array([0.8, 0.8])\n",
    "def FB0(P0):\n",
    "    xi0, c0 = P0\n",
    "    xb0 = XB_Default_B(mu0,c0)\n",
    "    F1 = ((taue - taui) * c0/r - I) * beta1 + (beta1 - 1) * (1 - taue) * xi0/(r - mu0) + (beta1 - beta2) * ((taui - taue) * c0/r - alpha * (1 - taue) * xb0/(r - mu0)) * (xi0/xb0)**beta2\n",
    "    F2 = (1 - taui) * c0/r + ((1 - alpha) * (1 - taue) * xb0/(r - mu0) - (1 - taui) * c0/r) * (xi0/xb0)**beta2 - p\n",
    "    return [F1, F2]\n",
    "xi0, c0 = fsolve(FB0, y0)\n",
    "xb0 = XB_Default_B(mu0,c0)\n",
    "E0 = E_Equity_Mu(xi0,mu0,c0)\n",
    "D0 = D_Debt_Mu(xi0,mu0,c0)\n",
    "phi0 = Phi_GCost_Mu(xi0,mu0,c0)\n",
    "F0x = F_Option_Mu(xi0,mu0,c0,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SEB(P0):\n",
    "    xis, cs = P0\n",
    "    Es = E_Equity_Mu(xis,mu0,cs) \n",
    "    F0s = F_Option_Mu(xi0,mu0,c0,xis) \n",
    "    Ps = P_FGap_Mu(xis,mu1,cs)\n",
    "    F1=Es-(I-Ps)-F0s \n",
    "    F2=Ps-p   \n",
    "    return [F1, F2]\n",
    "y0 = np.array([1, 1])\n",
    "ximax, cs = fsolve(SEB, y0) \n",
    "xis = min(ximax,xi1) \n",
    "xbs = XB_Default_B(mu1,cs)\n",
    "Es = E_Equity_Mu(xis,mu1,cs)\n",
    "Ds = D_Debt_Mu(xis,mu1,cs)\n",
    "phis = Phi_GCost_Mu(xis,mu1,cs)\n",
    "Fs = F_Option_Mu(xis,mu1,cs,x)\n",
    "Cs=(Fg-Fs)/Fg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:14.967697Z",
     "start_time": "2023-11-04T11:41:14.539517Z"
    }
   },
   "outputs": [],
   "source": [
    "tb0, tb1 = 0, 1.3  \n",
    "xlim0, xlim1 = tb0, 1 \n",
    "yb0, yb1 = xis, xi0 \n",
    "if kappa < 0:   \n",
    "    tb1=1.07 \n",
    "yb0=0.0000000001   \n",
    "M = 1000   #time \n",
    "N = 100    #space\n",
    "shape = [M, N] \n",
    "t, y = np.linspace(tb0, tb1, shape[0]), np.linspace(yb0, yb1, shape[1])\n",
    "dt, dy = t[1] - t[0], y[1] - y[0]\n",
    "T, Y = np.meshgrid(t, y, indexing='ij')\n",
    "Fpt = np.zeros(shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:15.936421Z",
     "start_time": "2023-11-04T11:41:15.496217Z"
    }
   },
   "outputs": [],
   "source": [
    "it1=np.where(xlim1-0.000001<T[:,0])\n",
    "it1=it1[0][0]\n",
    "nt=M-it1\n",
    "nd = 1\n",
    "ip=np.where(1.69999<Y[0,:])\n",
    "ip=ip[0][0]\n",
    "xp = int(Y[0,ip]*100)/100\n",
    "ip=np.where(1.5409999<Y[0,:])\n",
    "ip=ip[0][0]\n",
    "xp = round(Y[0,ip]*100)/100\n",
    "ipx1=np.where(x<Y[0,:])\n",
    "ipx1=ipx1[0][0]\n",
    "xpx1 = round(Y[0,ipx1]*100)/100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:18.331710Z",
     "start_time": "2023-11-04T11:41:18.042768Z"
    }
   },
   "outputs": [],
   "source": [
    "pixt = 1 / (1 + (1 - pi0) * ((x / Y) ** (psi / sigma ** 2)) * np.exp(kappa * T) / pi0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:22.164362Z",
     "start_time": "2023-11-04T11:41:18.334704Z"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure() \n",
    "ax = fig.add_subplot(projection='3d')\n",
    "ax.plot_surface(T[:-nt,1:-nd], Y[:-nt,1:-nd], pixt[:-nt,1:-nd], cmap='viridis')\n",
    "ax.set(xlabel='Elapsed time t', ylabel='Cash flow x')  \n",
    "plt.title(f'Posteriori probability $\\pi(x,t)$')\n",
    "plt.show()\n",
    "fig_name = 'pi_3D_surface_sigma' + str(int(sigma*100))\n",
    "fig.savefig(f'{fig_name}.eps', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if sigma < 0.2:  \n",
    "    pixt0 = pixt\n",
    "    ip0 = ip\n",
    "    T0 = T\n",
    "    nt0=nt\n",
    "if sigma > 0.2: \n",
    "    pixt1 = pixt\n",
    "    ip1 = ip\n",
    "    T1 = T\n",
    "    nt1=nt\n",
    "if np.isclose(sigma,0.2): \n",
    "    pixt2 = pixt\n",
    "    ip2 = ip\n",
    "    T2 = T\n",
    "    nt2=nt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'pixt0' in dir() and 'pixt1' in dir() and pixt0.shape==pixt1.shape:\n",
    "    nyt = 4\n",
    "    fig = plt.figure()\n",
    "    plt.subplot(2, 1, 1)\n",
    "    plt.plot(T0[:-nt0,ip0], pixt0[:-nt0,ip0], '-b', lw=3)\n",
    "    ylim0, ylim1 = plt.ylim() \n",
    "    plt.ylim([ylim0, ylim1])\n",
    "    ytk0=np.arange(ylim0, ylim1, step=(ylim1-ylim0)/nyt)\n",
    "    ytk0=[0.5940, 0.5945, 0.5950,  0.5955, 0.5960]\n",
    "    plt.yticks(ytk0)\n",
    "    plt.legend(['$\\sigma=0.15$'])\n",
    "    plt.title(f'$\\pi(x,t)$ with $x$ = {xp}')\n",
    "    plt.grid(True, 'both')\n",
    "    plt.subplot(2, 1, 2)\n",
    "    plt.plot(T1[:-nt1,ip1], pixt1[:-nt1,ip1], '--r', lw=3)\n",
    "    ylim0, ylim1 = plt.ylim() \n",
    "    plt.ylim([ylim0, ylim1])\n",
    "    ytk1=np.arange(ylim0, ylim1, step=(ylim1-ylim0)/nyt)\n",
    "    ytk1 = [0.5346, 0.5349, 0.5352, 0.5355, 0.5358]\n",
    "    plt.yticks(ytk1)\n",
    "    plt.legend(['$\\sigma=0.25$'])\n",
    "    plt.xlabel('Elapsed time t')\n",
    "    plt.grid(True, 'both')\n",
    "    plt.show()\n",
    "    fig_name = 'pi_t2'\n",
    "    fig.savefig(f'{fig_name}.eps', bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:42.303053Z",
     "start_time": "2023-11-04T11:41:23.301808Z"
    }
   },
   "outputs": [],
   "source": [
    "VM = np.zeros_like(Y)  \n",
    "VM_Strike = np.zeros_like(Y)\n",
    "Ep = np.zeros_like(Y)\n",
    "for i in range(shape[0]):   \n",
    "    tV = T[i]  \n",
    "    yV = Y[i]  \n",
    "    def CouponV_Pool(PV):\n",
    "        cpgV = PV        \n",
    "        ppiV = 1 / (1 + (1 - pi0) * ((x / yV) ** (psi / sigma ** 2)) * np.exp(kappa * tV) / pi0)\n",
    "        xbp0V = (beta2 / (beta2 - 1)) * ((r - mu0) / r) * cpgV\n",
    "        xbp1V = (beta4 / (beta4 - 1)) * ((r - mu1) / r) * cpgV\n",
    "        D0V = (1 - taui) * cpgV / r + ((1 - alpha) * (1 - taue) * xbp0V / (r - mu0) - (1 - taui) * cpgV / r) * (yV / xbp0V) ** beta2\n",
    "        D1V = (1 - taui) * cpgV / r + ((1 - alpha) * (1 - taue) * xbp1V / (r - mu1) - (1 - taui) * cpgV / r) * (yV / xbp1V) ** beta4\n",
    "        phipV = 1 - (r * ((1 - ppiV) * D0V + ppiV * D1V)) / ((1 - taui) * cpgV)        \n",
    "        FV = (1 - phipV) * (1 - taui) * cpgV / r - p\n",
    "        return FV\n",
    "    cV=np.ones(yV.shape)*0.8\n",
    "    cpgV = fsolve(CouponV_Pool, cV)  \n",
    "    ppiV = 1 / (1 + (1 - pi0) * ((x / yV) ** (psi / sigma ** 2)) * np.exp(kappa * tV) / pi0)\n",
    "    xbp0V = (beta2 / (beta2 - 1)) * ((r - mu0) / r) * cpgV\n",
    "    xbp1V = (beta4 / (beta4 - 1)) * ((r - mu1) / r) * cpgV\n",
    "    Ep1V = (1 - taue) * (yV / (r - mu1) - cpgV / r + (cpgV / r - xbp1V / (r - mu1)) * (yV / xbp1V) ** beta4)\n",
    "    Dp0V = (1 - taui) * cpgV / r + ((1 - alpha) * (1 - taue) * xbp0V / (r - mu0) - (1 - taui) * cpgV / r) * (yV / xbp0V) ** beta2\n",
    "    Dp1V = (1 - taui) * cpgV / r + ((1 - alpha) * (1 - taue) * xbp1V / (r - mu1) - (1 - taui) * cpgV / r) * (yV / xbp1V) ** beta4\n",
    "    phipV = 1 - (r * ((1 - ppiV) * Dp0V + ppiV * Dp1V)) / ((1 - taui) * cpgV)\n",
    "    VM[i] = np.maximum(Ep1V - phipV * (1 - taui) * cpgV / r - (I - (1 - taui) * cpgV / r), 0)\n",
    "    if np.isclose(yb0, 0):\n",
    "        VM[i,0] = 0\n",
    "    VM_Strike[i] = phipV * (1 - taui) * cpgV / r + (I - (1 - taui) * cpgV / r)\n",
    "    Ep[i] = Ep1V\n",
    "    if i==0: \n",
    "        if np.isclose(kappa,0) and i==0:   \n",
    "            VFpk0x = (Ep1V - phipV * (1 - taui) * cpgV / r - (I - (1 - taui) * cpgV / r)) * (x / yV) ** beta3\n",
    "            imax = np.argmax(VFpk0x)\n",
    "            Xik0 = yV[imax]\n",
    "            VFpk0V = (Ep1V - phipV * (1 - taui) * cpgV / r - (I - (1 - taui) * cpgV / r)) * (yV / Xik0) ** beta3\n",
    "            Phik0 = phipV[imax]\n",
    "            Fpk0 = VFpk0x[imax]\n",
    "            Ppik0 = ppiV[imax]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:42.446657Z",
     "start_time": "2023-11-04T11:41:42.306034Z"
    }
   },
   "outputs": [],
   "source": [
    "F0V = np.zeros_like(Y)\n",
    "F1V = np.zeros_like(Y)\n",
    "FsV = np.zeros_like(Y)\n",
    "VFpk0 = np.zeros_like(Y)  \n",
    "for i in range(shape[0]):   \n",
    "    yV = Y[i]  \n",
    "    FsV[i] = (Es - phis * (1 - taui) * cs / r - (I - (1 - taui) * cs / r)) * (yV / xis) ** beta3          \n",
    "    F0V[i] = (E0 - phi0 * (1 - taui) * c0/r - (I - (1 - taui) * c0/r)) * (yV / xi0)**beta1\n",
    "    F1V[i] = (E1 - phi1 * (1 - taui) * c1/r - (I - (1 - taui) * c1/r)) * (yV / xi1)**beta3\n",
    "    if np.isclose(kappa,0):\n",
    "        VFpk0[i] = VFpk0V        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:42.838654Z",
     "start_time": "2023-11-04T11:41:42.732471Z"
    }
   },
   "outputs": [],
   "source": [
    "aa = 0.5 * dt * (sigma**2 * np.arange(1, N - 1)**2 - (r - d) * np.arange(1, N - 1))\n",
    "bb = 1 - dt * (sigma**2 * np.arange(1, N - 1)**2 + r)\n",
    "cc = 0.5 * dt * (sigma**2 * np.arange(1, N - 1)**2 + (r - d) * np.arange(1, N - 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:42.949377Z",
     "start_time": "2023-11-04T11:41:42.841646Z"
    }
   },
   "outputs": [],
   "source": [
    "F1V = np.flipud(F1V)\n",
    "FsV = np.flipud(FsV)\n",
    "VM = np.flipud(VM)\n",
    "F0V = np.flipud(F0V)\n",
    "v = np.zeros((M,N))\n",
    "v[0,:] = VM[0,:]  \n",
    "v[1:,0] = FsV[1:,0] \n",
    "v[1:,-1] = Ep[1:,-1]*np.exp(-d * np.arange(1, M) * dt) - VM_Strike[1:,-1] * np.exp(-r * np.arange(1, M) * dt) \n",
    "for i in range(1, M):  \n",
    "    v[i,1:N-1] = bb * v[i-1,1:N-1] + cc * v[i-1,2:N] + aa * v[i-1,0:N-2]\n",
    "def Cvam(vam0):\n",
    "    vam2 = np.empty_like(vam0)\n",
    "    if np.isclose(kappa, 0): \n",
    "        vam2[0,:] = VFpk0[0,:]  \n",
    "    else:  # Give an initial guess\n",
    "        if kappa > 0: \n",
    "            vam2[0,:] = (FsV[0,:]  + VM[0,:] )/2.03 \n",
    "        elif kappa < 0: \n",
    "            vam2[0,:] = (F1V[0,:]  + VM[0,:] )/2.1 \n",
    "    vam2[1:,0] = FsV[1:,0] \n",
    "    vam2[1:,-1] = Ep[1:,-1]*np.exp(-d * np.arange(1, M) * dt) - VM_Strike[1:,-1] * np.exp(-r * np.arange(1, M) * dt) \n",
    "    for i in range(1, M):  \n",
    "        vpayoff = VM[i,1:N-1]\n",
    "        vam1 = bb * vam2[i-1,1:N-1] + cc * vam2[i-1,2:N] + aa * vam2[i-1,0:N-2]\n",
    "        vam2[i,1:N-1] = np.maximum(vam1, vpayoff) \n",
    "    return vam2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:43.058748Z",
     "start_time": "2023-11-04T11:41:42.952944Z"
    }
   },
   "outputs": [],
   "source": [
    "v_init = np.zeros((M,N))\n",
    "vam = successive_approx(Cvam, v_init)\n",
    "v = np.flipud(v)\n",
    "vam = np.flipud(vam)\n",
    "F1V = np.flipud(F1V)\n",
    "FsV = np.flipud(FsV)\n",
    "VM = np.flipud(VM)\n",
    "F0V = np.flipud(F0V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:43.431404Z",
     "start_time": "2023-11-04T11:41:43.199558Z"
    }
   },
   "outputs": [],
   "source": [
    "Ytmx = np.zeros(shape[0])\n",
    "Yimx = np.zeros(shape[0],dtype=np.int8)\n",
    "Fpt = np.zeros(shape[0])\n",
    "for i in range(0, it1+1): \n",
    "    vpayoff = VM[i,1:N-1]\n",
    "    ## handle the points where payoff and vam are close to 0, considering numerical errors\n",
    "    ig0 = 0\n",
    "    ig0 = np.nonzero(  Ep[i,1:N-1] > VM_Strike[i,1:N-1] )[0] \n",
    "    ig0 = ig0[0]\n",
    "    diff_gf = VM[i,ig0+1:N-1]-vam[i,ig0+1:N-1]\n",
    "    diff_gf0i=np.nonzero(  np.isclose(diff_gf, 0)  )\n",
    "    imax = diff_gf0i[0] \n",
    "    if len(imax) == 0: \n",
    "        imax0 = N-1\n",
    "    else:\n",
    "        if kappa > 0:   \n",
    "            imax0 = imax[0] \n",
    "        if kappa < 0:   \n",
    "            # touching point inside, considering numerical errors\n",
    "            imax0 = np.int32(imax[0] + (imax[-1] - imax[0] + 1)*2/3) \n",
    "        imax0 = imax0 + ig0+1\n",
    "    Ytmx[i] = yV[imax0]\n",
    "    Yimx[i] = imax0\n",
    "    Fpt[i] = vam[i,imax0]\n",
    "Xipt = Ytmx\n",
    "# considering numerical errors\n",
    "_ , uni = np.unique(Xipt[:-nt], return_index=True) \n",
    "uni = np.sort(uni)\n",
    "if uni.size==1:\n",
    "    uni = range(M-nt)\n",
    "u1t = np.ones(np.size(uni))  \n",
    "if kappa < 0:   \n",
    "    uni = np.hstack((uni, M-nt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:49.057111Z",
     "start_time": "2023-11-04T11:41:48.113681Z"
    }
   },
   "outputs": [],
   "source": [
    "u=vam\n",
    "Phipt = np.zeros(shape[0])\n",
    "Ppit = np.zeros(shape[0])\n",
    "Cpt = np.zeros(shape[0])\n",
    "Cst = np.zeros(shape[0])\n",
    "Cptx1 = np.zeros(shape[0])\n",
    "Cstx1 = np.zeros(shape[0])\n",
    "Coupt = np.zeros(shape[0])\n",
    "Fst = np.zeros(shape[0])\n",
    "F0t = np.zeros(shape[0])\n",
    "F1t = np.zeros(shape[0])\n",
    "Ximint = np.zeros(shape[0])\n",
    "Xip0t = np.zeros(shape[0])\n",
    "for i in range(shape[0]):   \n",
    "    tt = T[i,0]   \n",
    "    yt = Xipt[i]\n",
    "    yi = Yimx[i]\n",
    "    def Coupont_Pool(P0):\n",
    "        cpgt = P0        \n",
    "        ppit = Pit_2Val(yt,tt,x)\n",
    "        Pt = P_FGap_Pit(yt,ppit,cpgt)\n",
    "        Ft = Pt - p\n",
    "        return Ft\n",
    "    c0=0.8\n",
    "    cpgt = fsolve(Coupont_Pool, c0)  \n",
    "    if np.isscalar(cpgt) == False:\n",
    "        cpgt = cpgt[0] \n",
    "    ppit = Pit_2Val(yt,tt,x)\n",
    "    phipt = Phi_GCost_Pit(yt,ppit,cpgt)\n",
    "    if np.isclose(kappa,0) == False:\n",
    "        def xip0t_equ(P0):\n",
    "            xip0, cp0 = P0\n",
    "            Ep0 = E_Equity_Mu(xip0,mu0,cp0) \n",
    "            F0p0 = F_Option_Mu(xi0,mu0,cp0,xip0) \n",
    "            ppit = Pit_2Val(xip0,tt,x)\n",
    "            Pp0 = P_FGap_Pit(xip0,ppit,cp0)\n",
    "            F1=Ep0-(I-Pp0)-F0p0  \n",
    "            F2=Pp0-p  \n",
    "            return [F1, F2]\n",
    "        P0 = np.array([1, 1])\n",
    "        xip0t, cp0t = fsolve(xip0t_equ, P0)  \n",
    "        Xip0t[i] = xip0t\n",
    "\n",
    "        def ximint_equ(P0):\n",
    "            xim, cm = P0\n",
    "            E1m = E_Equity_Mu(xim,mu1,cm) \n",
    "            Pm = P_FGap_Mu(xim,mu1,cm)\n",
    "            ipxim=np.where(xim<Y[i,:]) \n",
    "            ipxim=ipxim[0][0] - 1 \n",
    "            uOption = u[i,ipxim]\n",
    "            F1=E1m-(I-Pm)- uOption  \n",
    "            F2=Pm-p  \n",
    "            return [F1, F2]\n",
    "        if kappa > 0:\n",
    "            P0 = np.array([1.2, 1]) \n",
    "        else:\n",
    "            P0 = np.array([xip0t, 1]) \n",
    "        ximint, cmt = fsolve(ximint_equ, P0) \n",
    "        Ximint[i] = ximint\n",
    "        if kappa < 0 and i>0: \n",
    "            # little num errors on the right end\n",
    "            if Ximint[i] < Ximint[i-1]:\n",
    "                Ximint[i] = Ximint[i-1]\n",
    "\n",
    "    Fst[i] = FsV[i,yi]      \n",
    "    F0t[i] = F0V[i,yi] \n",
    "    F1t[i] = F1V[i,yi]\n",
    "    Cst[i]=(F1t[i]-Fst[i])/F1t[i]\n",
    "    Cpt[i]=(F1t[i]-Fpt[i])/F1t[i] \n",
    "    Cstx1[i]=(F1V[i,ipx1]-FsV[i,ipx1])/F1V[i,ipx1] \n",
    "    Cptx1[i]=(F1V[i,ipx1]-u[i,ipx1])/F1V[i,ipx1] \n",
    "    Ppit[i] = ppit\n",
    "    Phipt[i] = phipt\n",
    "    Coupt[i] = cpgt\n",
    "\n",
    "    if np.isclose(kappa,0) and i==int(shape[0]*2/3):   \n",
    "        Xitk0 = Xipt[i]\n",
    "        Phitk0 = Phipt[i]\n",
    "        Fptk0 = Fpt[i]\n",
    "        Ppitk0 = Ppit[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if kappa > 0:\n",
    "    locg = 'lower center'\n",
    "else: \n",
    "    locg = 'lower right'\n",
    "lw=3\n",
    "fig = plt.figure()\n",
    "plt.hlines(xi0, xlim0, xlim1, 'c', '--', lw=lw)\n",
    "plt.plot(T[uni,0], Xipt[uni], '-b', lw=lw) \n",
    "plt.hlines(xi1, xlim0, xlim1, 'k', ':', lw=lw)\n",
    "plt.plot(T[uni,0], Xip0t[uni], '-.r', lw=lw) \n",
    "plt.hlines(xis, xlim0, xlim1, 'g', '--', lw=lw)\n",
    "# consider numerical errors\n",
    "x_sm, y_sm = smxy(T[uni,0], Ximint[uni]) \n",
    "plt.plot(x_sm, y_sm, '-m', lw=lw) \n",
    "ylim0, ylim1 = plt.ylim() \n",
    "plt.ylim([ylim0, ylim1]) \n",
    "if kappa < 0:\n",
    "    ic=np.where(xis<y_sm) \n",
    "    ic=ic[0][0]\n",
    "    plt.plot(x_sm[ic], y_sm[ic], 'sm', lw=lw)\n",
    "    plt.vlines(x_sm[ic], ylim0, y_sm[ic], 'm', ':', lw=lw)\n",
    "    plt.legend(['FB with low $\\\\bar{x}_i(\\mu_L)$', 'PE with learning $x_{ip}(t)$', 'FB with high $\\\\bar{x}_i(\\mu_H)$', 'min PE with learning $x_{ip0}(t)$', '(max) SE with high $\\\\hat{x}_{i}=x_{imax}$', 'min SE with high $x_{imin}(t)$'], fontsize=8, loc=locg)\n",
    "plt.xlabel('Elapsed time t')\n",
    "plt.ylabel('Investment threshold')\n",
    "plt.title('Equilibrium on $(x,t)$-Plane, ' + '$\\sigma$ = ' + str(sigma))\n",
    "plt.show()\n",
    "fig_name = 'Xi_dyn' + '_sigma' + str(int(sigma*100))\n",
    "fig.savefig(f'{fig_name}.eps', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:49.072068Z",
     "start_time": "2023-11-04T11:41:49.060101Z"
    }
   },
   "outputs": [],
   "source": [
    "if np.isclose(kappa,0) == False:\n",
    "    if kappa > 0:  \n",
    "        Xipt0 = Xipt\n",
    "        Cpt0 = Cpt\n",
    "        Cptx10 = Cptx1\n",
    "        Ppit0 = Ppit\n",
    "        Phipt0 = Phipt\n",
    "        Fpt0 = Fpt\n",
    "        Coupt0 = Coupt\n",
    "        uni0 = uni\n",
    "        T0 = T\n",
    "        ip0 = ip\n",
    "        vamp0 = vam\n",
    "        Ximint0 = Ximint\n",
    "        Xip0t0 = Xip0t\n",
    "        nt0=nt\n",
    "    if kappa < 0: \n",
    "        Xipt1 = Xipt\n",
    "        Cpt1 = Cpt\n",
    "        Cptx11 = Cptx1\n",
    "        Ppit1 = Ppit\n",
    "        Phipt1 = Phipt\n",
    "        Fpt1 = Fpt\n",
    "        Coupt1 = Coupt\n",
    "        uni1 = uni\n",
    "        T1 = T\n",
    "        ip1 = ip\n",
    "        vamp1 = vam\n",
    "        Ximint1 = Ximint\n",
    "        Xip0t1 = Xip0t\n",
    "        nt1=nt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures - 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T11:41:55.353418Z",
     "start_time": "2023-11-04T11:41:53.036234Z"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure() \n",
    "ax = fig.add_subplot(projection='3d')\n",
    "ax.plot_surface(T[:-nt,1:-nd], Y[:-nt,1:-nd], u[:-nt,1:-nd], cmap='viridis')\n",
    "ax.set(xlabel='Elapsed time t', ylabel='Cash flow x') \n",
    "plt.title(f'Option value \\n $\\hat F_p(x,t;\\mu_H)$')\n",
    "plt.show()\n",
    "fig_name = 'vam_3D_surface_sigma' + str(int(sigma*100))\n",
    "fig.savefig(f'{fig_name}.eps', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figures - 2 sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'Xipt0' in dir() and 'Xipt1' in dir() and np.isclose(kappa,0)==False and Xipt0.shape==Xipt1.shape:\n",
    "    fig = plt.figure()\n",
    "    plt.plot(T0[:-nt0,ip0], vamp0[:-nt0,ip0], '-b', T1[:-nt1,ip1], vamp1[:-nt1,ip1], '--r', lw=3)\n",
    "    ylim0, ylim1 = plt.ylim() \n",
    "    plt.ylim([ylim0, ylim1])\n",
    "    # plt.xlim([xlim0, xlim1])\n",
    "    plt.yticks(np.arange(ylim0, ylim1, step=(ylim1-ylim0)/15))\n",
    "    plt.legend(['$\\sigma=0.15$', '$\\sigma=0.25$'])\n",
    "    plt.xlabel('Elapsed time t')\n",
    "    plt.ylabel('$\\hat F_p(x,t;\\mu_H)$')\n",
    "    plt.title(f'Option value $\\hat F_p(x,t;\\mu_H)$ with $x$ = {xp}')\n",
    "    plt.grid(True, 'both')\n",
    "    plt.show()\n",
    "    fig_name = 'vam2_time'\n",
    "    fig.savefig(f'{fig_name}.eps', bbox_inches='tight')    \n",
    "    \n",
    "    fig = plt.figure()\n",
    "    plt.plot(T0[uni0,0], Xipt0[uni0], '-b', T1[uni1,0], Xipt1[uni1], '--r', lw=3)\n",
    "    plt.legend(['$\\sigma=0.15$', '$\\sigma=0.25$'])\n",
    "    plt.xlabel('Elapsed time t')\n",
    "    plt.ylabel('Investment threshold')\n",
    "    plt.title('PE with learning')\n",
    "    plt.grid(True, 'both')\n",
    "    plt.show()\n",
    "    fig_name = 'Xipt2'\n",
    "    fig.savefig(f'{fig_name}.eps', bbox_inches='tight')\n",
    "\n",
    "    fig = plt.figure()\n",
    "    plt.plot(T0[uni0,0], Phipt0[uni0], '-b', T1[uni1,0], Phipt1[uni1], '--r', lw=3)\n",
    "    plt.legend(['$\\sigma=0.15$', '$\\sigma=0.25$'])\n",
    "    plt.xlabel('Elapsed time t')\n",
    "    plt.ylabel('Guarantee cost')\n",
    "    plt.title('PE with learning')\n",
    "    plt.grid(True, 'both')\n",
    "    plt.show()\n",
    "    fig_name = 'Phipt2'\n",
    "    fig.savefig(f'{fig_name}.eps', bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "445px",
    "left": "161px",
    "top": "135.688px",
    "width": "179px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "346.844px",
    "left": "951px",
    "right": "20px",
    "top": "120px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
